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ABSTRACT 

We present the results of a numerical study comparing photometric and physical 
properties of simulated z = 6 — 9 galaxies to the observations taken by the WFC3 
instrument aboard the Hubble Space Telescope. Using cosmological hydrodynamical 
simulations we find good agreement with observations in color-color space at all stud- 
ied redshifts. We also find good agreement between observations and our Schechter 
luminosity function fit in the observable range, M uv < —18, provided that a moderate 
dust extinction effect exists for massive galaxies. However beyond what currently can 
be observed, simulations predict a very large number of low-mass galaxies and evolving 
steep faint-end slopes from oil = —2.15 at z — 6 to olt = —2.64 at z = 9, with a de- 
pendence of oc (1 + z) 59 . During the same epoch, the normalization <fi* increases 
and the characteristic magnitude M* v becomes moderately brighter with decreasing 
rcdshift. We find similar trends for galaxy stellar mass function with evolving low- mass 
end slope from an = —2.26 at z = 6 to ajvf = —2.87 at z = 9, with a dependence 
of |qm| oc (1 + z) 65 . Together with our recent result on the high escape fraction of 
ionizing photons for low-mass galaxies, our results suggest that the low-mass galaxies 
are important contributor of ionizing photons for the reionisation of the Universe at 
z > 6. 

Key words: method : numerical — galaxies : evolution — galaxies : formation - 
cosmology : theory 



1 INTRODUCTION 

When and how galaxies formed throughout the history of 
the Universe is one of the most fundamental questions of 
astronomy and astrophysics. As technology improves, as- 
tronomers are able to push the redshift frontier of galaxy 
observation and attempt to obtain a glimpse of galaxy 
formation processes at high redshift. In particular, up- 
grades to the Hubble Space Telescope (HST) via the addi- 
tion of the Wide Field Camera 3 (WFC3) have expanded 
our view of the Universe to reach redshifts z > 7, and 
dozens of new high-redshift candidates have been identi- 
fied ( e.g.. lYan et al l |2009l ; lOesch et al.l 120091; lOuchi et all 
120091; iBouwens et alj|20ld : IWilkins et al.ll201lh . For exam- 
ple. iBouwens et al.l ( 20ld ) identified 113 z ~ 7 — 8 Lyman- 
Break candidates by color-color selection with AB magni- 
tudes M uv — 26 — 29 and faint-end luminosity slopes of 
a L = -1.94±0.24 at z ~ 7 and a L = -2.00±0.33 at z ~ 8. 
While these new observational results are very excit- 



ing, candidate galaxies at z > 7 are very faint and it is 
still difficult to obtain high-resolution spectra to estimate 
their physical properties such as stellar masses and star for- 
mation rates (SFR). Therefore it would be useful to ob- 
tain predictions from theoretical models on their spectro- 
photometric properties. In particular, cosmological hydro- 
dynamic simulations based on the standard concordance 
A cold dark matter (CDM) model can directly simulate 
the dynamical evolution of dark matter and baryons from 
the initial condition that is consistent with the observa- 
tions of cosmic microwave background. With the imple- 
mentation of radiative heating and cooling, one can also 
follow the thermal properties of intergalactic medium and 
collapse of gas clouds into galaxies. Much work has been 
done for c omparing simulated galaxies at z < 6 to observa- 
tions (e.g.jNapmine et al.ll2004l. l2005allbl: [Night et al.ll2006l ; 
iFinlator et alj 2006 ; iRobertson et al.ll2007 f). but the compar- 
isons of photometric properties of galaxies at z > 7 between 
simulations and observations have not been performed ade- 
quately yet. 
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One of the important missing information from theo- 
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Run 


Box size 


N p 


TUDM 




e 


-2" end 




(h-iMpc) 


(Dark Matter, Gas) 


(h-^M Q ) 




(ft-ikpc) 




N400L10 


10.0 


2x400 3 


9.37x10 s 


1.91x10 s 


1.0 


2.75 


N400L34 


33.75 


2x400 3 


3.60X10 7 


7.34x10 s 


3.38 


1 


N400L100 


100.0 


2x400 3 


9.37x10 s 


1.91x10 s 


6.45 





N600L100 


100.0 


2x600 3 


2.78x10 s 


5.65X10 7 


4.30 






Table 1. Simulation parameters used in this paper. The parameter N p is the number of gas and dark matter particles; moM and m gas 
are the particle masses of dark matter and gas; e is the comoving gravitational softening length; z cn d is the ending redshift of each 
simulation. 



retical predictions is the precise determinations of evolu- 
tion of galaxy luminosity function (LF) and stellar mass 
function (GSMF) at z > 6. Some of the earlier numerical 
work have suggested very steep faint-end slope of ai ^ — 2 
at z « 6 (|Night et all l200fj ; iFinlator et aUbood l. but its 
evolution has not been quantified yet at z > 6 using cos- 
mological hydrodynamic simulations. This is an important 
topic, because the number density of early dwarf galaxies has 
significant implications for the reionisation of the Universe 
at z > 6. It is thought that small dwarf galaxies, that are 
currently beyond detection, may have playe d an important 
role in the reionisation of the Universe (e.g. , iMufioz fc Loebl 
l20ld : lTrenti et al.ll2010l ; lYauma et alj|2010h . 

In this paper, we compare the results of our cosmo- 
logical simulations with the recent WFC3 observations. In 
particular, we will determine the evolution of faint-end (low- 
mass end) slope oll (oim ) at z — 6 — 9 more precisely in our 
simulations. The paper is organized as follows: Section[2]will 
briefly describe the simulation and setup that was used to 
obtain our data; Section [3] will outline the methods used 
in processing the photometric data; Sections [4][6] will detail 
assembled data from simulations; Section will discuss the 
implications of these findings to the epoch of reionisation. 
Finally in Section [9] we will summarise our findings. 



2 SIMULATIONS 

We use a modified version of the smoothed particle hydro- 
dynamics (SP H) code GADGET-3 (originally described in 
ISpringefeOOa ) . This m odified code includes radi ative cooling 
by H, He, and metals (|Choi fc Nagamindl2009l). heating by 
a unif orm UV bac kground of a mo d ified lHaardt fc Madaul 
(1 19961) spectrum (iKatz et all Il99rj ; iDave et all 1 19991 ). an 
lEisenstein fc Hul (l999h initial po wer spectrum, star for- 
mation via the "Pressure mod el" (Schave fc Dalla Vecchial 
|200S| ; IChoi fc Nagamind l20ld) . supe rnova feedback, sub- 
resolu tion model of multiphase ISM (Springel fc Hernquistl 
2003), and th e multicomponent varia ble velocity (MVV) 
wind model (|Choi fc Nagaminel 1201 ll ). We also explore 
the effect of a new UV background prescription by 
iFaucher-Giguere et al.l (2009) in Section 18.21 Our current 
simulations do not include the AGN feedback. 

We have shown earlier that the metal line cooling en- 
hances star formation acr oss all redshifts by about 10 — 30% 
jChoi fc Nagam inc 2009), and that the Pressure SF model 
suppresses star formation at high-redshift due to a higher 
threshold density for SF JChoi fc Nagamin el |2010|) with re - 
spect to the earlier model by Springel fc Hernquistl (2003). 



IChoi fc Nagaminel l|201ll) also showed that the MVV wind 
model, which is based on the momentum-driven wind, 
makes the faint-end slope of GSMF slightly shallower com- 
pared to the constant ve locity galactic wind model of 
ISpringel fc Hernquistl (|2003l ). One of the points of this paper 
is to report that, even with these improvements to the super- 
nova feedback models, our simulations continue to predict 
very steep faint-end slopes. 

It is know n that the collisional ionizat ion equilibrium 
(CIE) table of ISutherland fc Dopital (l993) overestimates 
the cooling rates compared to the case w hen the photoion- 
isation by the UVB i s take n into account jEfstathiou|[l992l ; 
IWiersma et"afl 120091 ). and IChoi fc Nagaminel (20091 ) have 
used this CIE table. We have compared our LF results in 
the two runs with and without the metal cooling, and found 
that there is little difference in the faint-end slope of the 
two simulations. Therefore this issue is not a problem for 
the present work. 

For this paper, simulations are setup with either 2 x 400 3 
or 2 x 600 3 particles for both gas and dark matter. Multi- 
ple runs were made with comoving box sizes of 10/i _1 Mpc, 
34/i _1 Mpc and 100/i _1 Mpc. These runs will be referred to as 
N400L10, N400L34 and N600L100, respectively. Simulation 
parameters are summarised in Table [1] The adopted cosmo- 
logical parameters are based on the most recent WMAP 
data: « m = 0.26, Oa = 0.74, O b = .044, h = 0.72, 
a 8 = 0.80, n s = 0.96 (Komatsu et alj|20ich . 

At each time step of the simulation, the gas particles 
that exceed the SF threshold density (nfjf = 0.6 cm -3 ) is 
allowed to spawn star particles with the probabili t y con - 
sistent with the intended SFR. I Choi fc Nagamind (2010T ) 
have shown that the Pressure SF model with the above 
nfh produces fa vorable results compared to the observed 
iKennicuttl (|l998l ) law, in particular at the low column den- 
sity end. Galaxies are then identified by a group-finder algo- 
rithm using a simplifi ed variant of the SUBFIND algorithm 
|Springel et al.ll200ll ). Collections of star and gas particles 
are grouped based on the baryonic density field. Properties 
such as stellar mass, SFR, metallicity and position, among 
others are recorded in galaxy property files. A more de- 
ta iled description of the group-finder process can be found 
in lNagamine et all (|2004h . 



3 METHOD 

Once we identify the individual galaxies, we calculate the 
spectral energy distributi on (SEP) of each galaxy by ap- 
plying the 2007 version of iBruzual fc Chariot! (2003) stellar 
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Figure 1. Panel (a) shows the observed- frame relative flux /„ of a simulated galaxy in the N400L34 run, redshifted from z = 6 to 
2 = 9 with no extincti on. The w edge-like feature on the blue-side is caused by the iMadaul jl995h prescription of IGM absorption. Panel 
(b) shows the effect of Calzctti (1997) extinction on the simulated galaxy spectrum at z = 7 for E(B — V) = 0.0,0.15, & 0.30. 



population synthesis model that includes the contribution of 
TP-AGB phase. We have repeated t h e calc ulation using the 
2003 version of iBruzual fc CharlotJ (|2003h . and confirmed 
that the results hardly change. This is expected because we 
are focusing our analyses on the rest-frame UV wavelengths 
range at z > 6, and the effect of TP-AGB appears only in 
the infra-red range. In the present work, nebular emission 
lines are not included in the SEDs, but our results are not 
affected because this work is mainly concerned about the 
comparisons in the rest-frame UV, where the optical emis- 
sion lines do not come into the filters. We will investigate the 
impact of nebular emission lines on the estimates of stellar 
masses of high-redshift galaxies in the subsequent work. 

After we calculate the intrinsic SEDs, we then ac- 
count for the dust extinction effect using the ICalzettil 
l|l997l ) extinction law. We adopt uniform constant values 
of E(B - V) = 0.0,0.075,0.10,0.15 & 0.30. The above 
extinction values are select ed based on the var ied data 
found in recent publications. iBouwens et~ai1 l|2009h claimed 
a trend of decreasing extinction approaching E(B — V) — 
at z ~ 6 based on the U V continuum slope. However 
Schacrcr fc de Barrosl (|2010h performed SED fits including 
the effect of nebular emission lines, and showed that dust at- 
tenuation of up to E(B - V) ~ 0.32 (A v ~ 1 for R v = 3.1) 
can be found for galaxies at z = 6 — 8, and that E(B — V) 
might be greater for higher mass galaxies. Our final choice 
for E(B — V) at each redshift is based on a visual match of 
the LF to the observed data range, staying within the above 
mentioned range. In Section [5] we also examine the effect of 
variable extinction models. 

The final step to achieving a spectrum as it would ap- 
pear to HST, is to factor in the absorption by the inter- 
galactic medi um (IGM). To do this we apply the model of 
lMadaul(|l995l ). Figure [TJi shows the relative flux /„ of a rep- 



resentative galaxy redshifted to z — 7, 8 & 9. The HST filter 
functions used by the observers are also shown on the X- 
axis: 1775, z 8 50, iios, J125, H 160 centered at 0.775, 0.850, 
1.05, 1.25, 1.60 fim, respectively, are some of the HST fil- 
ters strategically positioned to identify the Lyman break of 
high r edshift galaxies. Figure [Tp shows the effect of lCalzettll 
(|1997l ) extinction on the simulated galaxy spectrum at z = 7 
for E(B -V) = 0.0, 0.15, & 0.30. 

The observer-frame spectrum is then used to calculate 
the apparent magnitude of each galaxy for each HST filter, 
as well as the rest-frame UV magnitude centered at 1700A. 
This process gives us a rest-frame magnitude in each fil- 
ter that can then be assembled into color-color plots and 
luminosity functions (LFs) for direct comparison to those 
produced by observers. 



4 COLOR-COLOR DIAGRAMS 

In Figures(2]|3]&|4l we present color-color diagrams at z = 6, 
7 & 8, respectively. For each redshift we make direct com- 
parisons to the observed data. In each figure, the results 
from three simulations (N400L10, N400L34, N600L100) are 
shown separately, and the simulated galaxies are represented 
by blue, green, and red points for E(B — V) = 0.0, 0.15, & 
0.30, respectively. Galaxy evolution tracks are shown by the 
corresponding solid blue, red and green lines, which were 
produced by redshifting the average spectrum of simulated 
galaxies. By doing this we are able to show how the Lyman- 
break of the galaxy spectrum is passing through the HST 
filter set (Figure QJ,), transitioning into and out of the drop- 
out candidate selection area which has been defined by ob- 
servers. For example in Figure(2] 1775 — 2:850 for z = 6 galaxies 
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Figure 2. Color-color diagram comparing observed z ~ 6 candidates with our simulations of differing box sizes. The red, green and 
blue data points represent simulated galaxies with E(B — V) = 0.0, 0.15, &0. 30, respectively. For each extinction, the color track of 
a representative gala xy is shown with a solid line from 2 = 5 through z = 7 (out of plot r ange). The grey shad ed area, defined by 
( BouwensetalJ[2006), indicates the expected position of z ~ 6 objects. Black filled squares Bouwens et alll2006h and open triangles 
l Yan &; Win dhorst 2004) indicate observed z ~ 6 candidates. 




Figure 3. Same as Figure |21 but for 2^7. The grey shaded area, defined by lOesch et al ] j2009l ). indicates the expected position of 
2 ~ 7 objects. For each extinctio n, the color track o f a representative ga laxy is shown w ith a solid line from 2 = 5 through 2 = 8 (out of 
plot range). Black filled squares llOesch et al 1 l2009l) and open triangles llYan et alj|2009l) indicate observed z ~ 7 candidates. 




-0.5 0.5 1 1.5 -0.5 0.5 1 1.5 -0.5 0.5 1 1.5 



125 160 1S5 ISO 125 160 

Figure 4. Same as Figure [2] but for z ~ 8. The grey shaded area, defined by l|Bouwens et alj|2010h . indicates the expected position of 
2 ~ 8 objects. For each extinctio n, the color track of a representative gala xy is shown with a solid line from 2 = 6 through 2 = 9 (out of 
plot range). Black filled squares llBouwens et alj|2010l) and open triangles l l Van e t al. 2009) indicate observed z ~ 8 candidates. 
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Figure 5. This figure explores the origin of scatter in the color-color diagrams. Panel (a) shows the effect of varying dust extinction. 
Inset plot shows a gaussian distribution of E(B-V) around a center value of 0.15 with a a = 0.038. Panel (b) shows the effect of varying 
redshift. Inset plot shows a gaussian redshift distribution around z = 6 with a a = 0.25. Panel (c) shows the effect of varying dust 
extinction, redshift and adding a a = 0.125 gaussian distribution of photometric errors to the AB magnitudes in each filter. Inset plot 
shows distribution of AB magnitudes error. 



are substantially redder than at for z — 5, creating a thresh- 
old above which candidates are selected. 

In Figure [2] observed c andidate data points are repre- 
sented by solid b lack squares (IBouwens et alj2006h and open 
black triangles (|Yan fc Windhorst 20041 ). The grey shaded 
area indicates t he co lor-color selection criteria defined by 
iBouwens et al.l (|2006t ). Due to the abundance of observed 
points at this epoch, error bars were omitted for clarity. 

Figure [3] sh ows color-color plot s for z ~ 7 with observed 
data point s fromlOesch et al.l l|2009l , black filled squares) and 
lYan et all l|2009l . open triangles) . The grey shaded a rea for 
color-color selection is taken fro mlOesch et alj (120091) . Error 
bars were included only for the lOesch et al. | i|2009h data for 
clarity. Again we see good agreement with observed data. 

Figure [4] shows the c olor-color diagram for z ~ 8 with 
observed dat a points fro m Bouwens et al.1 (|2010l . black filled 
squares) and lYan et al.l (|2009l . open triangles). Expected po- 
sition for z ~ 8 obje cts is defined by the grey shaded area 
(IBouwens et alj|20ich . Error bars were included only for the 
IBouwens et al. ( 201ot ) data for clarity. Here too simulated 
galaxies show good agreement with observations. 

In all of Figures [2] [3] & [4] massive galaxies are on the 
lower left side of the data distribution for the simulated 
galaxies, because higher mass galaxies have higher SFR in 
our simulations and hence bluer. 

In Figure [S] we examine the tight grouping of galaxies 
in our color-color plots. Figure^ shows the effect of adding 
a ±0.15 scatter with a a — 0.038 gaussian distribution cen- 
tered around E(B — V) = 0.15. This addition does little to 
increase the scatter of the galaxies in color-color space. In 
Figure [5p we demonstrate the effect of uncertainties in pho- 
tometric redshift estimation by adding a ±0.50 scatter with 
a (7 = 0.25 gaussian distribution centered on z = 6. This ad- 
dition produces a vertical feature which runs throughout the 
selection box with some horizontal scatter at its top. Finally 
in Figure [5}; we add a ±0.20 scatter with a a = 0.125 gaus- 
sian distribution to the calculated UV magnitudes as well 
as including the previous scatter in E(B — V) and redshift. 
In this case, we reproduce a similar degree of scatter to the 
observed data. Our results indicate that the current obser- 



vational color-color selection region can select galaxies with 
a wide range of extinction values of E(B — V) = 0.0 — 0.30, 
and that the scatter in the data points are primarily domi- 
nated by the redshift scatter and photometric errors rather 
than the extinction variance. The tight grouping of galaxies 
seen in each of the color-color plots is because the simulated 
galaxies are extracted at exactly z = 6, 7, 8 and have no 
errors included in magnitude calculations. 



5 LUMINOSITY FUNCTIONS 

The rest-frame UV LFs of simulated galaxies at z = 6 — 9 
are shown in Figure [5] To represent a larger range of mag- 
nitudes, we construct composite LFs using three simula- 
tions with different box sizes, and connect the results at 
the points where they overlap. At the brightest end, the 
data correspond to the last bin in which we have data in 
the N600L100 run. At the faintest end, the minimum M uv 
corresponds to the minimum M s allowed by the mass res- 
olution in the N400L10 run (see Section [6}. This also cor- 
responds to the point at which the LF begins to turn over, 
indicating the resolution limit. For example, at z = 6, the 
magnitude coverage of each run are approximately from 
M uv « -23 to -20 mag (N600L100), -20 to -17.5 mag 
(N400L34), and -17.5 to -15mag (N400L10). These val- 
ues vary slightly with redshift. In all panels, three sets of 
results for E(B - V) = 0.0,0.10,0.30 at z = 6 & 7, and 
E(B - V) = 0.0,0.075,0.30 at z = 8 & 9 are shown. See 
Section [3] regarding these extinction choices. 

To explore the im pact of a varyin g extin ction, we utilize 
a model presented in iFinlator et all (120061 ): E(B - V) = 
9.0 x Z°' 9 +SE, where Z is the galaxy stellar metallicity and 
5E is a Gaussian scatter. This relationship is based on the 
mean trend between metallicity and E(B— V) found in SDSS 
at z ~ 0. In Figure[6ja) the result of this relationship applied 
to our simulated galaxies is shown by the cyan solid line. 
This relationship gives higher extinction for more massive 
galaxies, therefore our data is consistent with observations 
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Figure 6. Composite LFs at z = 6 - 9 using the N400L10 (green), N400L34 (blue), N400L100 (red) runs, as well as multiple ex- 
tinction values of E(B — V) = 0.0 (solid lines) , 0,075 (dotted line s), 0.10 (dashed lines), and 0.30 (dot-dashed). The solid cyan line 
represents a variable model for dust extinction IIFinlat or et al.||2006|). For 2 = 6. the yellow shade represents a range of observatio ns by 
l|Yan fc Windhorsdliool iBouwens et al.ll2004! : iDickinson et al.ll2004l : lMalhotra et al.ll2005l;|Bunker et al.ll2004l: iBeckwith et al.ll2006h. For 
z = 7, the yellow shaded area represents a range of obse rvations by jB*mrw ens**et*ldJl200^: lMcTure^r*alIl2'oiOl ; IOescn*^m*H20*loT " 20091 ; 
iMannucci et al]|2007i. IRichard et alj|2006h assembled in llBouwens et al.l|20ldV For 2 = 8, observed range represented in yellow is taken 
from llBouwens et al.ll20ld ) and includes data from jMcLure et al.ll2010l ; ICastellano et al.ll2010l ; IBouwens et al.ll20ld) and there are no 
observed data points for 2 = 9. 



only at the very brightest end of the LF with this extinction 
model. 

The disagreement at fainter magnitudes between the 
cyan line and the yellow observed region has several poten- 
tial explanations. The Finlator's model for E(B — V) is based 
on the observed relationship at z ~ 0, and thus it may not 
apply at high redshifts. This model also utilizes metallic- 
ity taken from star particles, which is set to the same value 
as the parent gas particle at the time of formation. There- 
fore the metallicity of star particles in our simulation does 
not evolve with redshift. We have tested this model utiliz- 
ing the gas particle metallicity, and found that this change 
results in lower extinction values, however these values may 
be artificially lower since this mean metallicity is obtained 
by averaging over the gas particles including those at the 
outskirts of galaxies. Finally there is the possibility that if 
high redshift galax ies are mostly dust free as s uggested by re- 
cent observations ()Bouwens et al ] |2009l . l201ll ). then our sim- 
ulations are systematically overproducing galaxies by ~ 0.5 
dex. This could be true if molecular gas cooling plays an 



important role in early star formation (see Section 18.21 for 
further discussion). However the consistent overproduction 
of galaxies is not a feature found in our GSMF, as we will 
discuss further in Section [6] and Figure flOl 

At z = 6 (top left panel of Figure |6l) , the yellow shade 
repres en ts a range of o bserva tio ns by lYan fc Windhorst! 
J2004I); IBouwens et all J2004); IDickinson et all J2004l; 



iMalhotra et alj |2005l ); lBunker et al.l (|2004T i; IBeckwith et all 
(2006). At this redshift we show good agreement with ob- 
served data. 

At z = 7, the yellow shaded a rea re pre sents a range of 
observations by IRichard et al.l J2006I); Mannucci et al.l 

Oesch et ail l|2009t ): 



<f2007h: 



Bouwens et alJ 



ll200Sl): 



iMcLure et al.l (120101 ): lOesch et all (|2010T l. Here too wc 
see good agreement with observations for the N400L34 and 



N600L100 run. 



At z 



8, the yellow shaded area re presents a range 



a rep: 

of obs e rvations bv IMcLure et aD ()2010h : ICastellano et al] 
l|20irj) ; IBouwens et alj (|2010l ). Again our simulation shows 
good agreement with observations. 



© 0000 RAS, MNRAS 000, 000-000 



Steep Faint-end Slopes of Galaxy Mass and Luminosity Functions at z > 6 7 




Figure 7. Best-fit Schechter LFs at z = 6 — 9 for the simulated galaxies. The simulation LFs (same as Figure |6j are shown by data 
points together wi th Poisson error bars . The yellow shaded areas for observed data are the same as in Figure |5] Schechter best-fits to 
the observed data llBouwens et alj|20ict dashed lines) is shown together with the 1-sigma error range (dotted grey lines). 



At z — 9, we only show the results from the N400L10 
and N400L34 runs as there are no objects found at this 
redshift in the N600L100 run due to low resolution. Also 
there has been no published data for robust observed LF 
estimates at z = 9. 

5.1 Schechter Function Fit 

We perform x 2 fit analysis of the ISchechterl (|l976l ) func- 
tion to the composite simulated galaxy LFs. The Schechter 
function, <^(Af uv ) = dn/dM uv , as a function of magnitude is 
given by 

<HM UV ) = 0l(O.41nlO)lO°' 4(M -- M - )(1+a ^ 

exp(-10°- 4 [ M --Muv] ); (1) 

where <f>* , M* v and ql are the normalization, characteristic 
magnitude, and faint-end slope. 
We minimize the expression 

X^E^t (2) 

by systematically adjusting each of the three Schechter pa- 
rameters, where Oi is our simulated data set, Ei is the cal- 
culated value of the above Schechter function, and JVbins = 



z a L log(0*) M* v x 2 



6 


-2.151 


24 
.15 


-3.461 


.29 
.37 


-21.151 


.53 

.53 


0.04 


7 


-2.30l 


.28 
.18 


-3.741 


.32 
.41 


-20.821 


.61 
.56 


0.06 


8 


-2.51+ 


.27 
.17 


-4.30+. 


.32 
.412 


-21. oo! 


.59 
.50 


0.11 


9 


-2.64+. 


.12 
.13 













Table 2. Best-fit Schechter LF parameters at z = 6 — 9 for the 
simulated galaxies. Some parameters at z = 9 are unavailable due 
to lack of data at the bright end from the N600L100 run. 



16, 14, 15, & 10 for z = 6, 7, 8, & 9, respectively. For z = 9 a 
simple power-law was used, because the N600L100 run did 
not contain any galaxies at this redshift due to limited reso- 
lution. The 1-sigma error of the parameters are determined 
by fixing two of the best -fit parameter s and varying the third 
from Xmin to Xmin + 1 (| Andraell20T0l ) . 

Figure [7] shows our x 2 best-fits to the Schechter func- 
tion along with the o bservational fits (long-dashed line) by 
iBouwens et al.l (120101) . The dotted l ine re presents the 1- 
sigma error range of bouwens et all l|2010l ) fit. The yellow 
shaded areas are the same as in Figure [6] For the Schechter 
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fits, we used the simulation data sets with E(B — V) = 0.10 
for z = 6 & 7, and E(B - V) = 0.075 for z = 8 & 9, as they 
agree with the yellow shade best. The best-fit Schechter pa- 
rameters are given in Table [2] 

Both our Schechter fit and the simulation LFs agree well 
with the yellow observed range at the bright end. However 
our best-fit Schechter LFs are noticeably steeper than the 
observed one at the faint end. This is caused by numerous 
faint galaxies with Muv > — 18 mag in our simulation, which 
lie beyond the current detection limit. Our results suggest 
that future surveys with increased sensitivity (such as the 
JWST) will detect a large number of faint galaxies. Although 
at different redshifts, this trend of finding more galaxies be- 
low the previous magnitude limit seems to be occur ring in 
recent observations. For example. iDrorv et a l. (2009) found 
an upturn toward a steep faint-end slope of ql ~ 1.7 at lower 
mass at z < 1, when they probed do wn to fainter limi t s than 
before in the COSMOS field. Also, iReddv fc Steidell (|2009T ) 
found a steeper slope of ql = —1.74 at z ~ 2 — 3 than before 
using a larger sample of > 2000 spectroscopic redshifts and 
~ 31000 LBGs. 

We summarise the evolution of Schechter parameters 
of simulated galaxies in Figure [S] The yell ow shades show 
the ra nge of observed data points taken from lBouwens et al.l 
(2010). In the top left panel, our simulation shows a clear 
evolution of ql from ctL = —2.15 at z — 6 to oil = —2.51 
at z = 8, and our slopes are in general steeper than the 
current observations by Aql — 0.3. At z = 9, we used a 
power-law fit to the incomplete data set (missing bright-end 
in the N600L100 run), therefore the z — 8 & 9 data points 
are connected by a dashed line. 

Our simulations suggest a steepening trend of 

\a L \ = 0.68(1 + z) - 59 , or (3) 
a L = -0.17(z - 6) - 2.16. (4) 

This steepening of «i t owards higher r e dshift is in congru- 
ence with the result of iBouwens et all (|201ll ). which finds 
a similar redshift evolution. As we will show in the next 
section, this steepening of ql is driven by the simultane- 
ous steepening of galaxy stellar mass function. In a ACDM 
universe, we expect a larger number of lower mass galaxies 
at higher redshifts, and these low mass galaxies gradually 
merge to form more massive galaxies. Therefore the steep- 
ening of a with increasing redshift is naturally expected in 
our simulations. 

The characteristic magnitude M* v in our simulation 
does not show much evolution, except that it brightens by 
about 0.3 mag from z — 7 to z = 6. The simulation M* v 
is brighter than the observed range by 0.5 — 1 mag. The 
normalization parameter <f) L of simulated galaxies increases 
significantly by about an order of magnitude from z = 8 to 
z — 6 with a redshift dependence of 

^(x(l + Z )-°- 42 . (5) 

The simulation <j> L is lower than the observed one by about 
0.5 dex. 

One of the possible reasons for our lower <j> L and 
brighter M^ v than the observational estimates is that our 
LF does not have a distinct exponential break at the bright- 
end, and the simulation LF seems more like a single power- 
law with a slight curvature at the bright-end. Therefore we 
are somewhat forcing a Schechter fit, which could push the 



Mu V to brighter side cf>* L to lower side. Future simulations 
with AGN feedback may make the bright-end break more 
pronounced, and the values of above two parameters might 
approach the observed range. 



6 GALAXY STELLAR MASS FUNCTIONS 

Galaxy stellar mass functions (GSMFs) in our simulations 
at z = 6 — 9 are shown in Figure [9] As with the luminosity 
function, data from three simulations have been combined 
to show a continuous mass function. The yellow shaded 
area shows the obse rved mass functions for z — 6 & 7 
l|Gonzalez et alj|20ich , which w as derived direc t ly fro m the 
luminosity function obtained bv lBouwens et all (|2010l ). The 
lower mass limit of M s = 10 6 ' 8 M Q for the N400L10 run 
was determined by converting the faint-end luminosity cut- 
off of Muv = —15 to M a via the relationship between M uv 
and log (M s ) in our simulations (Figure [10)) . In order to stay 
within the mass resolution limits of each run, we take the 
following mass range from each run: M s = 10 6 - 8 - 1O 8 M 
(N400L10), 1O 8 -1O 9 M (N400L34), > W 9 M Q (N600L100). 
This ensures that, based on the star particle mass for each 
run (Table [T|, we will have a sufficient number of particles 
in each galaxies. For example, the N400L10 run has a star 
particle mass of 9.55 x lO 4 Af0, which is half of the gas par- 
ticle mass. Therefore a galaxy of 1O 6 ' 8 M0 would contain 
approximately 66 star particles. Simulated galaxies would 
also contain gas and dark matter particles, so the actual 
particle number of simulated galaxies would be significantly 
larger. Further discussion of the simulation resolution limits 
can be found in Section [8. II 

The Schechter function <j>(M s ) = ln(10)M s $(M a ) = 
dn/dlog M s as a function of galaxy stellar mass M s is de- 
fined as 

W ) = 0alO)A(g) exp(-^), (6) 

where tj)^, Ms , and om are the normalization, characteristic 
stellar mass, and low-mass (faint-end) slope, respectively. 
We perform the same \ 2 analysis for these three Schechter 
parameters. The best-fit parameters are given in Table [3] 
Similarly to the LF, we find that the GSMF steepens from 
um = —2.26 at z — 6 to um = —2.57 at z = 8, with an 
evolutionary trend of 

\a M \ = 0.63(1 + z)°- 65 , or (7) 
a M = -0.20(z - 6) - 2.26. (8) 

Figure [5] includes the results of GSMF Schechter fits, 
showing a clear trend that galaxies become more massive 
and numerous from z = 9 to z = 6 as they grow with time. 
The characteristic number density tj>* M and the associated 
M* increase with decreasing redshift, evolving as 

<&oc(i+*r°- 46 , (9) 

and 

M* oc (l + z)- ' 20 . (10) 

These trends are also summarised in the bottom right panel 
of Figure [5] clearly showing the evolution of GSMF from 
z = 9 to z — 6. 

Interestingly the slope om is steeper for the GSMFs 
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GSMF 




Figure 8. Evolution of Schechter parameters for both LF (blue) and GSMF (red) as functions of redshift. The dashed line connecting 
2 = 8 data points to z = 9 indicates that the 2 = 9 data set is incomplete at the bright end. 1-sigma error bars are shown for all data 
points, see 15.11 for detail of error analysis. Red data points are offset for clarity in a and <j>* plots. Dotted red and blue lines indicate 
(I + 2) 7 fits to the data in each plot with t he exception of M*,, wh ich shows only a week trend with redshift. Yellow shaded area represents 
the observed fit parameter range found in lBouwens et all d201ll) with the black dashed line indicating the linear fit to that data. 
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log(<?>*) 



log(M s *) x 2 
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25 
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-2.431 
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-2.571 


15 
17 


-4.691 


.32 
.34 
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0.18 
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-2.871 


12 
12 











Table 3. Best-fit Schechter parameters for simulated GSMFs at 
2 = 6 — 9. The evolution of these parameters are also summarised 
in Figure[8] Some parameters at 2 = 9 are unavailable due to lack 
of data at the massive end from the N600L100 run. 



than for the LFs by Aa ~ 0.1. This can be explained by 
examining the relationship in our simulations between M uv 
and log (M s ) (Figure [TO)). Through a linear fit to the median 
data point in each M uv bin, we find 



log(M,) = /3Muv + C, 



(11) 



where (3 « -0.35 and C = 1.90,1.75,1.59,1.58 at z = 
6, 7, 8, 9, respectively. Since the faint-end slopes of the 
Schechter luminosity and mass functions (Eqs. [T] and [6]) 



are —0.4(1 + lvl) and (1 + qm), respectively, the relation- 
ship between the two slopes can be derived from dn = 
(j)(M uv )dM uv = $(M s )dM s as follows: 



0(1 



-0.4(1 + Q L ) 



(12) 



Therefore plugging in /3 = —0.35 and solving for cxm results 
in 



a M = 0.14 + 1.14ol, 



(13) 



which clarifies the relationship between the two slopes. Uti- 
lizing this relationship, we are able to recover either am or 
ql to within < 0.5 sigma from one another. The exact value 
of f3 must depend on the details of the stellar IMF and star 
formation histories, and we plan to investigate the evolution 
of ft as a function of redshift in more detail in future work. 

We see in Figure [9| that our simul a tion r esults and the 
observational result of Gonzalez et al.l l|2010l ) do not agree 
well at both faint and bright end, with Gonzalez's GSMF be- 
ing shallower than our simulation result. This can be clearly 
explained through different M uv vs. log(M s ) relationship, as 
we show in Figure 1101 The Gonzalez's relationship (shown 
as blue long-dashed line; log M s = — 0.68M UV + const; their 
Fig. 1) is much steeper than what our simulation suggests, 



© 0000 RAS, MNRAS 000, 000-000 



10 Jaacks et al. 




Figure 9. Composite GSMFs at z = 6 — 9 with best-fit Schechter LFs. The best-fit Schechter pa rameters are summaris ed in Table [3] 
Yellow shaded regions at z = 6 and z = 7 indicate observational derived mass functions taken from Gonzalez et al.l (2010) . 



and for a given value of M uv , they infer a higher M s for a 
brighter M uv , and a lower M s for a fainter M uv . We note 
that they derived this relationship from z = 4 data set, 
and the limited data at z — 6 does not seem to follow 
this relationship tightly. This difference in the slope of M uv 
vs. log(Af s ) relationship is directly reflected in the differ- 
ence in the GSMF slope. If we translate M uv into luminos- 
ity, then our s i mulat ion predicts Luv oc Mf' 875 , whereas 
iGonzalez et all l|201Ch obtained Luv oc Af* -7 at z ~ 4. We 
plan to study the evolution of relationship between M uv , M s 
and SFR in our simulations in our subsequent work. 

In the next section, we discuss the implications of these 
steep mass/luminosity functions in our simulations on the 
reionisation of the Universe. This is an important subject 
that is closely related to our current work, as the source of 
ionising radiation is considered to be coming from these low 
mass, star forming galaxies at z > 6. 



7 IMPLICATIONS FOR REIONISATION & 
MINIMUM HALO MASS FOR HOSTING 
FAINT GALAXIES 

In order to quantify the contribution of the faint, low-mass 
galaxies to the reionisation of the Universe, we examine the 
SFR density (SFRD) from galaxies with different stellar 



masses at z > 6. Figure [TT] shows the SFRD as a func- 
tion of redshift broken down into galaxy stellar mass bins of 
M s = lO 6 ' 8 - 10 s M (red), 10 8 - 10 9 M Q (blue), > 1O 9 M 
(black), and the total SFRD (dashed dark green). We also 
show the observational esti mates taken f rom Bouwens et al.l 
(|2010l . black squares) and lOuchi et al.l (|2009l . magenta tri- 
angles). The fractional contribution from each mass range 
is listed in Table [4] and its redshift evolution is summarised 
in Figure [T2l 

The dotted black lines indicate the SFRD required to 
maintain reionisation (|Madau et al.|[l999h : 



2 x 10" 



C 

I esc 



l + z 
10 



(14) 



which is dependent on redshift and the ratio of IGM clump- 
ing factor (C) and escape fraction (/ esc ) of ionizing photons 
from galaxies. Although the exact value for C is unknown, a 
range of value s has been estimated fro m numeri cal simula- 
tions : C = 30 dGnedin fc Ostrikedl997h. C = 10 (jlliev et all 
120061) . and C = 3-6 jPawlik et all 1 200a). 

In our recent work, Yaiima et al.l 1 20101 ) estimated the 
values of f eac as a function of halo mass in our simula- 
tions, by performing ray-tracing radiative transfer calcula- 
tions of ionizing photons from star-forming galaxies. They 
found that f esc decreases with increasing halo mass. In or- 
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Figure 10. Relationship between rest-frame UV magnitude and 
galaxy stellar mass in our simulations indicated by the black con- 
tours (N600L100, N400L34 and N400L10 from upper left to lower 
right) for z = 6. Median values in each 0.5 mag bin are shown by 
the black and red triangles for 2 = 6 and 8, respectively, with 
a least-square linear fit to the median points (solid black and 
red short-dashed lines, with a slope of m —0.35). Points greater 
than our M uv limit of —15 are exclude from t he fit. The blue 
long-d ashed line is the observational result from [Gonzalez et al.l 
||2010|) with a slope of -0.68. 



der to obtain a representative value of / esc , we compute 
the average (/esc), weighted by the halo number density, us- 
ing [Sheth]&JIbrmen| (|2002h halo mass function, and obtain 
(fesc) = 0.42. Here we use this value to calculate C for each 
assumed value of C// esc . 

For comparison, the yellow and cyan shaded ar- 
eas in Figure 1111 repres ent a range of model predic- 
tions from iMunoz fc Loebl l|2010l ). They developed a semi- 
analytic model of star-forming galaxies using extended 
Press-Schechter theory, and examined the effect of minimum 
halo mass that can host a galaxy on the UV LF. They found 
that the best-fit UV LF to the current observations can be 
obtained with a minimum halo mass of Mh, m in ~ 1O 9 ' 4 M0. 
The cyan shade indicates the SFRD derived using this best- 
fit minimum halo mass and varying L1500/ SFR. The same 
is done considering Afh.mm = lO 8 A/0, shown by the yel- 
low shade. By comparison, our simulation is more consis- 
tent with the case of Mh, m i n ~ 10 s M , as confirmed by the 
direct data in our simulation (see the later discussion in Sec- 
tion [8]3] and Figure [16)) . In our N400L10 run which resolves 
the lowest halo mass among the three simulations, we find 
Afh,min ~ 10 8 ' 2 M©, corresponding to the virial temperature 
of T v i r » 10 4 K. This is expected, because the cooling curve 
used in our simulation cuts off at around this temperature. 

Using the co nditional luminosity function method, 
iTrenti et all \20ld) found a similar trend of a steepening 
faint-end slope from ai = —1.74 at 2 = 6 to oll — —1.99 at 



Figure 11. Star formation rate density (SFRD) as a function of 
redshift. The result from multiple galaxy stellar mass ranges are 
shown: 10 6 - 8 < M 3 /M Q < 10 8 (dotted red), 10 8 < M a /M e < 
10 9 (short dashed blue) and M s > 1O 9 M (dashed-dot black). 
The total sum of all mass ranges is shown as the dark green 
dashed line, and the fractional contribution from each mass range 
can be found i n Tab le [4] Observational es timates take n from 
iBouwens et al. I <201Cft (black squares) and lOuchi et all (|2009l) 
(magenta triangles). The yellow and cyan shaded areas represent s 
a range of model estimates presented in iMuhoz fc Loebf feOldl 
Dotted black contour lines represent minimum SFRD requ ired to 
keep the IGM ionized as defined by iMadau et all (|l99Sh using 
varying values for C and a constant value of f esc = 0.42. 



SFRD Fraction 



z 


10 6.8 _ 10 8 Mq 


10 8 - 1O 9 M 


> 1O 9 M 


6 


0.66 


0.25 


0.08 


7 


0.72 


0.24 


0.04 


8 


0.83 


0.15 


0.02 


9 


0.89 


0.11 




10 


1.00 







Table 4. Fractional contribution from each galaxy stellar mass 
range to the total SFRD found in Figure [TTI at z = 6 — 10. The 
redshift evolution of these values are summarised in Figure 1121 



z = 9, based on the evolution of dark matter halo mass func- 
tion. Given these steep faint-end slopes, they argued that 
galaxies between M uv — — 18 and M uv = — 10 contribute 
> 75% of the total luminosity density at z = 7 and produce 
enough ionizing photons to maintain reionisation assuming 
C/feac 25 (C — 5, / esc = 0.2). In our simulations, this 
fraction is even higher (> 90% for the same M uv range as 
above) due to our steeper ol- However, if we change the 
integration limit to — 18 < A/ uv < — 15 corresponding to our 
resolution limit, we find this fraction to be 74% at z = 7. 
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6 7 8 9 



z 

Figure 12. Fractional contribution from each galaxy mass range 
to the total SFRD found in Figure [TT] at z = 6 - 10. See Table g] 
for the numerical values. 



Z PL Ps 

-18 < Aiuv < -15 10 6 - 8 < M s < 1O 8 - 2 M 

6 0.64 0.71 

7 0.74 0.82 

8 0.83 0.89 



Table 5. Fractional contribution to the UV Luminosity density 
(Pl) and stellar mass density (p s ) by the faint-end galaxies be- 
low current detection limit, calculated by integrating the best-fit 
Schechter functions. Limits of integration correspond to the cur- 
rent observational limit (M uv = —18, M s = 1O 8 2 M0) and to the 
simulation resolution limit (M uv = —15, M a = 10 6 ' 8 Mq). See 
Sections [5] and [6] for details regarding resolution limits. 



See Table [5] for the fractional contribution by the faint-end 
galaxies to the total UV luminosity density and stellar mass 
density at z — 6 — 9. 

Examining the total SFRD found in our simulations 
(long-dashed line), we find that, depending on our assump- 
tions of the values for C, galaxies can maintain reionisation 
as early as z ~ 9.5 for C = 1, and as late as z ~ 6.5 
for C — 30. In our simulations, low-mass galaxies with 
M 3 = 10 6 ' 8 - 10 s M© (red line) and M uv > -18 are the 
primary contributor to the SFRD at z > 6 as shown in 
Table |3J and they are capable of achieving reionisation by 
z ~ 6. Note that our simulations do not include contribu- 
tions from Pop III stars. 



0.9 



0.8 - 



.0.7 - 



0.6 



0.5 




6.5 



7.5 



Figure 13. Fractional contribution to the total luminosity den- 
sity (pl) and stellar mass density (p a ) by the galaxies below the 
current detection limit of M uv = —18. The values are obtained 
by integrating the best-fit Schechter functions in the range of 
-18 < A/uv < -15 or M s = 10 6 ' 8 - 10 8 M©. See Table [5] for 
numerical values. 



8 DISCUSSIONS 

8.1 Numerical Resolution Issues 

In order to assess the resolution limit of the simulation, one 
may compare the particle mass and smoothing length to the 
Jeans mass and length. Assuming the SF threshold density 
(nfjf = 0.6 cm -3 ) and a temperature range of 500 — 3000 K, 
we find that the Jeans length and mass to be Aj ~ 300 — 
700 pc and Mj « 10 5 — 1O 7 M0. Compared to these numbers, 
the proper gravitational softening length of the N400L10 
run at z — 6 is e = 197 pc, which is shorter than Aj by 
100 — 500 pc. The gas particle mass of the same run (m gas = 
2.65 x 10 5 Af e with h = 0.72) is lower than Mj by a factor 

of w 1 - 100. 

However, iBate fc Burkertl (|l997l ) argued for a mass res- 
olution criteria that requires the local Jeans mass to be less 
than M res ~ 2iV no i g h, where -/V ne i g h is the number of parti- 
cles in the SPH kernel. For the N400L10 run, iVncigh = 33 
was used, resulting in M ICB ~ 6.6 x 1O 6 M0. This number is 
in-between the range of Mj given above, and it suggests 
that our N400L10 simulation cannot resolve the gravita- 
tional collapse of gas clouds properly at temperatures be- 
low ~2120 K, similarly to other cosmological SPH simula- 
tions of galaxy formation with similar box sizes and particle 
counts. This is one of the reasons why cosmological simu- 
lations need to adopt a sub-particle multiphase ISM model 
and form star particles out of a cold phase gas that co- 
exist with a hot phase gas with high effective temperature. 
Such sub -particle SF mode ls are designed to reproduce the 
observed Kennicuttl |l99ct) l aw in a statistical sense (e.g. , 
ISpringel fc Hernquistj 12003 ; ISchave fc Dalla Vecchial 120081 ; 
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Figure 14. Luminosity functions at z = 6 from the runs of 
the same box size (10h~ 1 Mpc) and varied resolutions: 2 X 144 3 
particles (short-dashed blue line), 2 X 216 3 particles (dash-dot red 
line) and 2 X 400 3 particles (long-dashed black line). The dotted 
black line represents the faint-end slope of a = —2.15. 



IChoi fc Nagaminell2010h , even though the simulations can- 
not resolve the collapse of the cold gas clouds. 

Despite of all the discussion above, one would have to 
perform a convergence study to be really confident that the 
results are unaffected by the resolution limit. To show that 
our results on the faint-end are not resolution dependent, 
we compare the LFs created from runs with seed particle 
counts of 2 x 400 3 (long dash black), 2 x 216 3 (dash-dot red) 
and 2 x 144 3 (short dashed blue) all in a 10/i _1 Mpc box 
size in Figure [T4l Here we see that the faint-end slope from 
M uv = —18 to M uv = — 15 remains largely unchanged in the 
N400L10 run. 

8.2 UVB Radiation and H 2 Cooling 

Among the results presented in this paper, the most striking 
one is the very steep faint-end slopes in both galaxy mass 
and luminosity functions at z > 6, and we find a greater 
number of low-mass galaxies than what the current obser- 
vations suggest. While it is possible that the current obser- 
vations are missing these faint galaxies at high-redshift and 
that the upcoming JWST will detect these faint sources, we 
should also consider the uncertainties and limitations in the 
treatment of star formation in our numerical simulations. 

In our simulations, stars are allowed to form once the 
gas density exceeds the th reshold density of nfjf = 0.6 cm~ 3 
(see lNagamine et al.l2O10l . for justification of this value), and 
star particles are generat e d acco rding to the SF law matched 
to the local iKennicuttl (|l998l ) law bas ed on the cooling 
funct ions with metal line contributions (|Choi fe Nagaminej 
2010). If the star formation efficiency decreases with increas- 
ing redshift in the real Universe, it is possible that our sim- 



Figure 15. Comparison of z = 7 LFs using different UVB pre- 
scriptions. The long-dashe d black line repr esents the N 400L10 run 
using a modified Haardt fc Madaul l ll996l) spectrum (Ka tz et al.l 
11999 iDave et all Il999l). The lo ng-short dashed red line is 
with iFaucher-Giguere et all (|2009| . FG-UVB) model in the same 
N400L10 run. The dash-dot blue line is a N400L10 run using 
FG-UVB and a variable star formation threshold density (VTH) 
which has a (1 + z) dependency to mimic the possible effect of 
molecular hydrogen cooling on star formation. The short-dashed 
black line (with a faint-end slope of a = —2.01) indicates the 
Schechter function fit to observed data (cy an shade) along with l a 
fit parameter errors (dotted black lines) by Bo uwens et al. N2011I) . 
Magenta triangles represent z = 7 LF with no dust extinction 
from the lKuhlen et al l ll201ll ) KMT09 simulation, which used an 
SF model based on the H2 mass. Filled magenta circles approxi- 
mate the results of applying a dust extinction of E(B — V) = 0.10. 



ulations are overproducing stars at high redshift. This could 
happen, for example, if star formation is controlled primar- 
ily by the amount of molecular gas, which is dependent on 
the gas metallicity. In the early Universe, ISM and IGM may 
not be enriched with metals and dust sufficiently yet, which 
leads to a lower molecular mass densities and lower SFRD 
i|Gnedin fc Kravtsov|[20lb1 ; iKuhlen et al.ll201lh . 

To investigate the possible impact of H2 cooling on our 
result, we experiment with a simple model varying star for- 
mation threshold (VTH), in which the value of nfjf is in- 
creased by a factor of (1 + z) with increasing redshift. The 
result of this run is shown by the dash-dot blue line in Fig- 
ure 1151 As the SF threshold density is increased and the 
star formation is limited to regions with higher and higher 
density, the number of low mass galaxies is decreased sig- 
nificantly at z = 7, showing a much flatter faint-end slope 
than our fiducial runs. 

In Figure [15] we also compare o ur fiducial and VTH 
model to the KMT09 simulation by IKuhlen et all (|201 ll ). 
which utilizes an SF model based on H2 mass (magenta tri- 
angles). Although covering only small part of the total LF, 
we can clearly see that without dust extinction correction, 
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log M h [M ] log M h [M ] 



Figure 16. Contour plot of absolute UV flux corresponding to 
AB magnitude at lOpc plotted against dark matter halo mass 
for each run N400L10 (left), N400L34 (middle) and N600L100 
(right). The luminosity for each halo is a summation of all galaxies 
found in the halo. Blue dashed line show a least square fit to the 
data. Solid re d line represents the results of a semi analytic model 
developed bv lTrenti et al. Yellow shaded region shows a 

range of predict ions based on semi-analytic calculations found in 
iLee et alj <200Sh . 



Figure 17. Composite dark matter halo mass function (HMF) 
at 2 = 6 (red), 2 = 7 (black) and 2 = 8 (green) assembled at each 
redshift from the N600L100, N400L34 and N400L10 runs. The 
dashed line at each redshift represents predicted halo ma s s func - 
tion based on theoretical work found in ISheth fc Tormenl d2002f> . 
Blue solid line indicates a slope= —2.00 offset for comparison. 



this model overpredicts the observed rest-frame UV LF at 
z — 7, similarly to our simulations. When we apply a dust 
extinction effect with E(B — V) = 0.10 approximately to 
the Kuhlen et al.'s data (a shift of ~1 mag similarly to our 
simulations), their result (filled magenta circles) agree with 
observations very well just like our result. This suggests that 
the simulations utilizing a SF model based on H2 mass pro- 
duce a very similar LF to that presented here, except the 
sharp turnover at M uv ~ — 18 seen for the Kuhlen et al.'s 
data. We are currently in the process of implementing a sim- 
ilar SF model based on H2 mass in our simulations, and we 
will report the results in a subsequent paper (Thompson, et 
al., in preparation). 



8.3 Origin of Steep Faint End Slope 

In order to compare our result with some of the semi- 
analytic models based on the halo occupation model, in Fig- 
ure [16] we examine the relationship between the rest-frame, 
monochromatic flux /1500 at 1500 A and the halo mass Mh 
for all the halos in our simulations at z = 6. The dashed 
blue line is a least-square fit to all the data points shown by 
the contour, the solid red line i s the result of semi-analytic 
modeling bv lTrenti et alj J2010[). and the yello w shade is an- 
other similar model result by Lee et all l|2009f ). Trenti et al. 
calibrated their model to the observed LF at z = 6, and 
Lee et al. to the observed LF and correlation functions at 
z — 4,5, 6, with both assuming single galaxy occupation per 



halo. Lee et al.'s model implements a mass and luminosity 
threshold for galaxies, thus the yellow shade does not cover 
the entire range of our simulation. 

We see good agreement between these model results 
at M h > 10 10 M Q . Below M h = 10 10 Af Q our relationship 
begins to deviate from Trenti's model, being higher by ~ 0.5 
dex. This result implies that, in our simulations the low mass 
halos with < 10 10 Mq are more efficient at producing stars 
than in the model of Trenti et al., making our faint-end slope 
steeper than theirs and Bouwens et al.'s. 

To further validate our simulation results, in Figure [T71 
we compare the dark matter halo mass function (HMF) 
in our simulations (da ta points) to the analytic model of 
ISheth fc Tormenl (|2002l , dashed line). The HMF was assem- 
bled from our three runs in the same manner as the LF and 
GSMF. Here we find very good agreement with the model, 
indicating that we are producing proper number of halos in 
the mass range of 10 s < M h < 1O 12 M . 

Since we are producing correct number of halos, the 
above results suggest that our steeper than observed faint- 
end slope can be attributed to the fact that the low mass 
halos are producing stars more efficiently, rather than too 
many galaxies in a particular Af uv bin. This would have 
the effect of steepening the faint-end slope, since the devia- 
tion occurs primarily in the N400L10 run which controls the 
faint-end of the LF. Further work must be done to explain 
the exact baryonic physics of this effect. 
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9 CONCLUSIONS 

Using cosmological SPH simulations, we examined the colors 
of high-redshift galaxies, quantified the evolution of luminos- 
ity and mass functions at z = 6 — 9, compared our results 
to the latest WFC3 observational results, and examined the 
implications for the reionisation of the Universe. Following 
are the main conclusions of the present work: 

• We find good agreement with observations in color-color 
space at z = 6, 7, and 8. Most of the simulated galaxies 
with E(B — V) — 0.0 — 0.30 are consistent with the color 
selection criteria with a relatively tight distribution. This 
suggests that the galaxies selected by the current color se- 
lection criteria could have a wide range of extinction values 
with E(B - V) = 0.0 - 0.30. We also find that the scatter 
in the observed data on the color-color plane is more domi- 
nated by the redshift scatter and photometric errors, rather 
than the variance in extinction. 

• The rest-frame UV LF of simulated galaxies agree well 
with the current WFC3 observations when we assume uni- 
form dust extinctions of E(B — V) = 0.10 for z — 6,7 and 
E(B — V) = 0.075 for z — 8. T hese extinction values are con - 
sistent with those obtained bv lSchaerer fc de Barrosl (|2010h . 
who performed SED fitting including nebular emission lines. 
See Section [5] for more details on the extinction treatment. 

• We performed least \ 2 fits to simulation LFs using the 
three parameter Schechter luminosity and mass functions. 
The best-fit Schechter parameters are summarised in Ta- 
bles [2] and [3] Each of the three Schechter parameters, (q_l, 
M* v , and 4>*l) or («m, M* , and 4>* M ), are evolving with red- 
shift. The faint-end slope is steeper than the current ob- 
servational estimates with a < —2, and becomes steeper 
with \a L \ oc (1 + z)°' 59 and \a M \ oc (1 + z) 0S5 ■ The char- 
acteristic mass M* decreases towards higher redshift with 
M* oc (1 + z) -0 ' 20 . The characteristic magnitude M* v does 
not evolve very much, but it does get dimmer from z — 6 
to 2 = 7 by about 0.3 magnitude. The normalization <f>* 
decreases towards higher redshift as expected in the hier- 
archical structure formation model, with ci/, oc (1 + z)~ 0A2 
and <jf M oc (l + z)~ ' 45 . 

• We decomposed the total SFRD into three separate 
contributions from different galaxy mass ranges of M s = 
10 6 ' 8 - 10 8 M o , 10 s -10 9 M B and > 10 9 M B , in a similar fash- 
ion to the work bv lChoi fc Naeaminel l|201ll ). We find that 
galaxies with M s = 10 6 ' 8 — 10 s Mq are the primary contrib- 
utor to the total SFRD at z > 6 as shown in Figs. 1111 1121 
and Table [H and therefore to the ionizing photon budget as 
well. Our simulation suggests that these faint galaxies can 
reionise the Universe by z = 6 as long as the clumping fac- 
tor is C < 30. This conclusion is consistent wi th that of our 
earlier direct radiative transfer calculations bv lYaiima et al.l 
l|2010h . 

• We find that the contribution from low-mass galaxies 
with M s = 10 6 ' 8 - 10 s M Q to the total SFRD (and hence the 
ionizing photon budget) decreases with decreasing redshift, 
from 89% at z — 9 to 66% at z — 6, which is still dominant. 
This decreasing trend can also be seen in the contribution 
from faint-end galaxies ( — 18 < M uv < —15) to the total 
UV luminosity density pL and stellar mass density p 3 (see 
Table [5] and Figure [13)) . Whereas the fractional contribution 



from galaxies with M s = 10 s - 1O 9 M and > 1O 9 M are 
increasing gradually during the same epoch with much lower 
fraction (see Tableland Figure [T2"j). The fractional contri- 
bution from galaxies with M s > 10 9 Mq is particularly small 
with less than 10% at all redshifts of z = 6 — 9. 

In summary, the steep faint-end slopes of high-redshift 
galaxies at z > 6 seems to be a generic prediction of ACDM 
model, i f one calibrates th e SF model to match the locally 
observed iKennicuttl (1 19981 ) law or to the lower redshift LF 
results. In addition to the results of our simulations, the two 
simple semianalytic models obtain similarly ste ep faint-end 
slopes (jMunoz fc Loebll2010l : iTrenti et~aD l2010h . as well as 
the more detai l ed sem ianalytic model of galaxy formation by 
ILo Faro et al.l |2009). The existence of these low- mass galax- 
ies below the current detection limit can be tested soon by 
the JWST, and we will continue to investigate the possi- 
bility of decreasing SF efficiency towards high-redshift by 
implementing more sophisticated physical processes (e.g., 
metallicity dependent SF model based on H2 masses) in our 
future simulations. Given the significant contribution by the 
faint-end galaxies to the overall budget of SFRD, ionizing 
photons, luminosity and stellar mass densities, it is criti- 
cal to determine the number density and the nature of this 
population more accurately. 
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